program calcP_kappa

        use mConstants
	use mObject
        use mCosmology
        use mLPS
        use mP_kappa
        use mPSFittingFunctions
	use mChi
        use mMisc

        implicit none

        real(D), parameter :: L_LL = 1e3_D, &
                              L_UL = 1e5_D

        integer, parameter :: NUMPTS = 50

        real(D) :: l(NUMPTS), factor, lightSmithP, darkSmithP
        integer :: i

        type(tUniverse) :: darkUniverse
        type(tLPS) :: darkLPS
        type(tChi) :: lightChi, darkChi
        type(tSmith) :: lightSmith, darkSmith

        call new(lightSmith)
        call new( lightChi, universe = universe(lightSmith) )

        call new(darkUniverse, w_0 = -0.8_D, w_a = -0.6) ! A new universe with dark energy.
        call new(darkLPS, universe = darkUniverse, normalization = 'CB') ! Linear power spectrum with dark energy and CMB normalization.
        call new(darkSmith, darkLPS)
        call new( darkChi, universe = darkUniverse )

        write( 6, '(a)' ) '# Smith et al. convergence power spectra:'
        write( 6, '(a)' ) '# k(h/Mpc) P_Smith P_Smith_DE'

        call geometricPartition(l, L_LL, L_UL)

        do i = 1, NUMPTS

                lightSmithP = P_kappa( lightSmith, lightChi, l(i) )
                darkSmithP = P_kappa( darkSmith, darkChi, l(i) )

                factor = l(i) ** 2 / (2.0_D * PI)

                write( 6, fmt = '(3e12.4)' ) l(i), &
                                             factor * lightSmithP, &
                                             factor * darkSmithP

        end do
        
end program calcP_kappa
